Planteo del Problema

El problema que vamos a tratar de resolver es predecir el salario de un jugador de la NBA para la temporada 2019-2020 en base a sus estadísticas de juego durante la temporada 2018-2019.

Las técnicas de regularización son útiles para trabajar con conjuntos con gran cantidad de variables, las cuales pueden introducir variabilidad en las estimaciones de los parámetros.

Conjuntos de datos

Los datos provienen de la página Basketball Reference y fueron previamente trabajados por nosotros para obtener el formato actual.

Las variables del set son:

diccionario = read_csv("diccionario_terminos.csv")
diccionario

En el glosario de Basketball Reference pueden encontrar una descripción más exhaustiva de cada una de estas métricas.

# Los datos de salario son para la temporada 2019-2020
nba = read_csv(" nba_player_stats_salary_2019_2020.csv") %>% 
  rename(salary = mean_salary_2019_2020) %>% 
  mutate(Pos = str_remove(string = Pos, pattern = "\\-.*")) %>% 
  mutate_all(~replace(., is.na(.), 0))

glimpse(nba)

## Rows: 395
## Columns: 50
## $ Player  <chr> "Jaylen Adams", "Steven Adams", "Bam Adebayo", "LaMarcus Aldr…
## $ Pos     <chr> "PG", "C", "C", "C", "SG", "C", "PF", "SF", "SF", "PF", "PF",…
## $ Age     <dbl> 22, 25, 21, 33, 23, 20, 28, 25, 25, 30, 24, 34, 21, 24, 33, 3…
## $ Tm      <chr> "ATL", "OKC", "MIA", "SAS", "UTA", "BRK", "POR", "ATL", "MEM"…
## $ G       <dbl> 34, 80, 82, 81, 38, 80, 81, 48, 43, 25, 72, 10, 67, 81, 69, 8…
## $ GS      <dbl> 1, 80, 28, 81, 2, 80, 81, 4, 40, 8, 72, 2, 6, 32, 69, 81, 70,…
## $ MP      <dbl> 428, 2669, 1913, 2687, 416, 2096, 2292, 463, 1281, 322, 2358,…
## $ FG      <dbl> 38, 481, 280, 684, 67, 335, 257, 64, 150, 21, 721, 49, 183, 1…
## $ FGA     <dbl> 110, 809, 486, 1319, 178, 568, 593, 157, 276, 69, 1247, 121, …
## $ `FG%`   <dbl> 0.345, 0.595, 0.576, 0.519, 0.376, 0.590, 0.433, 0.408, 0.543…
## $ `3P`    <dbl> 25, 0, 3, 10, 32, 6, 96, 24, 9, 9, 52, 21, 67, 81, 145, 131, …
## $ `3PA`   <dbl> 74, 2, 15, 42, 99, 45, 280, 77, 34, 40, 203, 64, 202, 217, 43…
## $ `3P%`   <dbl> 0.338, 0.000, 0.200, 0.238, 0.323, 0.133, 0.343, 0.312, 0.265…
## $ `2P`    <dbl> 13, 481, 277, 674, 35, 329, 161, 40, 141, 12, 669, 28, 116, 1…
## $ `2PA`   <dbl> 36, 807, 471, 1277, 79, 523, 313, 80, 242, 29, 1044, 57, 202,…
## $ `2P%`   <dbl> 0.361, 0.596, 0.588, 0.528, 0.443, 0.629, 0.514, 0.500, 0.583…
## $ `eFG%`  <dbl> 0.459, 0.595, 0.579, 0.522, 0.466, 0.595, 0.514, 0.484, 0.560…
## $ FT      <dbl> 7, 146, 166, 349, 45, 197, 150, 26, 37, 12, 500, 15, 36, 89, …
## $ FTA     <dbl> 9, 292, 226, 412, 60, 278, 173, 35, 64, 16, 686, 22, 62, 102,…
## $ `FT%`   <dbl> 0.778, 0.500, 0.735, 0.847, 0.750, 0.709, 0.867, 0.743, 0.578…
## $ ORB     <dbl> 11, 391, 165, 251, 3, 191, 112, 24, 48, 18, 159, 9, 58, 27, 5…
## $ DRB     <dbl> 49, 369, 432, 493, 20, 481, 498, 60, 203, 36, 739, 45, 139, 1…
## $ TRB     <dbl> 60, 760, 597, 744, 23, 672, 610, 84, 251, 54, 898, 54, 197, 2…
## $ AST     <dbl> 65, 124, 184, 194, 25, 110, 104, 23, 128, 19, 424, 5, 47, 269…
## $ STL     <dbl> 14, 117, 71, 43, 6, 43, 68, 22, 54, 4, 92, 4, 46, 65, 91, 52,…
## $ BLK     <dbl> 5, 76, 65, 107, 6, 120, 33, 13, 37, 1, 110, 7, 22, 4, 21, 4, …
## $ TOV     <dbl> 28, 135, 121, 144, 33, 103, 72, 23, 58, 14, 268, 8, 55, 63, 1…
## $ PF      <dbl> 45, 204, 203, 179, 47, 184, 143, 48, 112, 25, 232, 32, 140, 1…
## $ PTS     <dbl> 108, 1108, 729, 1727, 211, 873, 760, 178, 346, 63, 1994, 134,…
## $ PER     <dbl> 7.6, 18.5, 17.9, 22.9, 7.5, 18.5, 13.2, 11.2, 12.8, 4.6, 30.9…
## $ `TS%`   <dbl> 0.474, 0.591, 0.623, 0.576, 0.516, 0.632, 0.568, 0.516, 0.569…
## $ `3PAr`  <dbl> 0.673, 0.002, 0.031, 0.032, 0.556, 0.079, 0.472, 0.490, 0.123…
## $ FTr     <dbl> 0.082, 0.361, 0.465, 0.312, 0.337, 0.489, 0.292, 0.223, 0.232…
## $ `ORB%`  <dbl> 2.6, 14.7, 9.2, 10.3, 0.8, 9.6, 5.3, 5.3, 4.1, 6.1, 7.3, 3.3,…
## $ `DRB%`  <dbl> 12.3, 14.8, 24.0, 19.8, 5.1, 24.0, 22.6, 13.9, 18.1, 12.5, 30…
## $ `TRB%`  <dbl> 7.4, 14.7, 16.6, 15.1, 3.0, 16.8, 14.2, 9.5, 11.0, 9.2, 19.3,…
## $ `AST%`  <dbl> 19.8, 6.6, 14.2, 11.6, 8.9, 7.9, 6.0, 6.9, 15.0, 7.8, 30.3, 2…
## $ `STL%`  <dbl> 1.5, 2.0, 1.8, 0.8, 0.7, 1.0, 1.4, 2.2, 2.1, 0.6, 1.8, 0.7, 1…
## $ `BLK%`  <dbl> 1.0, 2.4, 3.0, 3.4, 1.1, 4.5, 1.2, 2.4, 2.7, 0.3, 3.9, 2.0, 1…
## $ `TOV%`  <dbl> 19.7, 12.6, 17.1, 8.8, 13.9, 13.0, 9.7, 11.8, 16.0, 15.5, 14.…
## $ `USG%`  <dbl> 13.5, 16.4, 15.8, 26.9, 24.4, 15.9, 13.7, 17.2, 12.6, 12.0, 3…
## $ OWS     <dbl> -0.1, 5.1, 3.4, 6.4, -0.4, 4.4, 3.0, 0.1, 0.8, -0.3, 8.9, 0.1…
## $ DWS     <dbl> 0.2, 4.0, 3.4, 2.9, 0.4, 3.3, 2.8, 0.4, 1.8, 0.0, 5.5, 0.2, 1…
## $ WS      <dbl> 0.1, 9.1, 6.8, 9.3, 0.0, 7.6, 5.8, 0.4, 2.7, -0.2, 14.4, 0.3,…
## $ `WS/48` <dbl> 0.011, 0.163, 0.171, 0.167, 0.002, 0.175, 0.121, 0.043, 0.100…
## $ OBPM    <dbl> -3.8, 0.7, -0.4, 2.4, -4.2, 0.2, 0.1, -2.4, -1.6, -4.5, 6.3, …
## $ DBPM    <dbl> -0.5, 0.4, 2.2, -0.6, -2.1, 1.4, 0.6, -0.3, 2.4, -1.9, 4.1, -…
## $ BPM     <dbl> -4.3, 1.1, 1.8, 1.8, -6.3, 1.6, 0.7, -2.7, 0.8, -6.4, 10.4, -…
## $ VORP    <dbl> -0.2, 2.1, 1.8, 2.6, -0.5, 1.9, 1.5, -0.1, 0.9, -0.4, 7.4, -0…
## $ salary  <dbl> 131678, 25842697, 3454080, 26000000, 2429400, 2376840, 925800…

Existen 395 observaciones y 50 variables en el dataset.

Análisis Exploratorio

Gráfico de la relación entre la posición y el salario

Veamos como es la distribución de salarios según la posición en el juego. Se agregan las etiquetas de los jugadores que cobran mayores sueldos.

top_players = c("James Harden", "Stephen Curry", "Blake Griffin", "Chris Paul", "LeBron James", "Klay Thompson", "Jimmy Butler", "Gordon Hayward", "Kyle Lowry")
ggplot(nba, aes(Pos, salary, fill=Pos)) +
  geom_boxplot() +
  geom_text(aes(label=ifelse((salary>30500000) & Player %in% top_players,as.character(Player),'')),hjust=1.1,vjust=0, size=3) +
  theme_bw() +
  labs(title= "Boxplots: salarios y posición de juego", x="Posición", y="Salario")

Observamos que la distribución varía un poco entre posiciones y hay varios jugadores que son outliers según el criterio de Tukey.

Correlograma

Realizamos un correlograma entre todas las variables cuantitativas.

nba %>% 
  select_if(is.numeric) %>% # selección variables numéricas
  ggcorr(., layout.exp = 5, hjust = 1, size = 3.5, nbreaks = 5, color = "grey50") + # graficamos correlacion pearson
  labs(title='Correlograma de variables cuantitativas')

Observamos que existen relaciones de diversa magnitud y signo entre todas las variables.

GGpairs (algunas variables)

Seleccionamos algunas variables y vemos sus relaciones usando ggpairs.

¿Cómo es la relación del salario con las demás variables?

nba %>% 
  select(salary, Age, PTS, GS, DRB, TRB, AST, BLK) %>% 
  ggpairs() + theme_bw()

Observamos que la correlación de todas estas variables con el salario es positiva pero de diversa magnitud.

También se ve que existen correlaciones muy fuertes entre ciertas variables. Por ejemplo, entre TRB (Total de Rebotes) y DRB (Rebotes Defensivos).

Modelo Lineal

Vamos a probar un modelo lineal que incluya todas las variables (excepto al jugador y equipo) y obtener las estimaciones de los parámetros junto a su p-valor e intervalo de confianza.

Coeficientes estimados

Creamos el modelo y accedemos a información de los coeficientes e intervalos usando tidy.

# Eliminamos jugador y equipo
nba = nba %>% select(-c(Player, Tm)) 
# Modelo lineal
modelo_lineal = nba %>% lm(formula = salary~., data = .)
#Coeficientes
lineal_coef= modelo_lineal %>% tidy(conf.int=TRUE)

Graficamos los coeficientes estimados.

lineal_coef %>% filter(!is.na(estimate)) %>% 
  ggplot(., aes(term, estimate))+
  geom_point(color = "forestgreen")+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), color = "forestgreen")+
  geom_hline(yintercept = 0, lty = 4, color = "black") +
  labs(title = "Coeficientes de la regresión lineal", x="", y="Estimación e Int. Confianza") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90))

Graficamos los p-valores de mayor a menor para evaluar la significatividad individual de los coeficientes estimados.

lineal_coef %>% filter(!is.na(estimate)) %>% 
  ggplot(., aes(reorder(term, -p.value), p.value, fill=p.value))+
  geom_bar(stat = 'identity', aes(fill=p.value))+
  geom_hline(yintercept = 0.05) +
  labs(title = "P-valor de los regresores", x="", y="P-valor") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90)) + 
  scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )

Notamos que:

Hay ciertos coeficientes estimados que presentan una gran variabilidad pero la escala de las variables puede ocultarnos la verdadera variabilidad de los estimadores.
Pocas variables tienen coeficientes significativos: PosPG, TOV%, G y Age.
Existen cuatro variables cuyo coeficiente estimado es NA: 2P, 2PA, PTS y TRB.

lineal_coef %>% filter(is.na(estimate))

Cuando una variable se puede expresar como una combinación lineal de otra, el modelo lineal de R devuelve los valores de los coeficientes estimados como NA.

Para evitar los problemas que puede introducir la escala, reescalamos las variables con el comando scale y repetimos el ajuste para estos nuevos datos.

# Reescalamos las variables numericas
nba_scaled = nba %>% mutate_at(vars(-Pos), scale)
# Nuevo modelo lineal 
modelo_lineal_scal = nba_scaled %>% lm(formula = salary~., data = .)
lineal_coef_scal = modelo_lineal_scal %>% tidy(conf.int=TRUE)

lineal_coef_scal %>% filter(!is.na(estimate)) %>% 
  ggplot(., aes(term, estimate))+
  geom_point()+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), color = "forestgreen")+
  geom_hline(yintercept = 0, lty = 4, color = "black") +
  labs(title = "Coeficientes de la regresion lineal", subtitle="Variables escaladas", x="", y="Estimacion e Int. Confianza") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90))

lineal_coef_scal %>% filter(!is.na(estimate)) %>%
  ggplot(., aes(reorder(term, -p.value), p.value, fill=p.value))+
  geom_bar(stat = 'identity', aes(fill=p.value))+
  geom_hline(yintercept = 0.05) +
  labs(title = "P-valor de los regresores",subtitle="Variables escaladas", x="", y="P-valor") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90)) + 
  scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )

Ahora observamos que:

Los coeficientes cambian de valor.
Las mismas cuatro variables tienen coeficientes significativos: PosPG, TOV%, G y Age.
Las mismas cuatro variables tienen coeficiente estimado NA: 2P, 2PA, PTS y TRB.

lineal_coef_scal %>% filter(is.na(estimate))

Evaluación de los modelos

Obtemos la evaluación de ambos modelos ¿Cómo esperan que sean los valores de diagnóstico?

modelos_lineales = list(lineal = modelo_lineal, lineal_escalado = modelo_lineal_scal)
map_dfr(.x = modelos_lineales, .f = glance, .id="modelo") %>% 
  select(modelo, r.squared, adj.r.squared, p.value)

La alta cantidad de variables y la existencia de una alta correlación entre varias de ellas ocasionan que los coeficientes estimados tengan alta varianza y que muchos de ellos no sean significativos en términos estadísticos. Las técnicas de regularización pueden ayudarnos a mejorar esta situación.

Regularización

La libreria glmnet nos permite trabajar con modelos ridge, lasso y elastic net. La función que vamos a utilizar es glmnet(). Es necesario que le pasemos un objeto matriz con los regresores y un vector con la variable a explicar (en este caso los salarios).

Fórmula:

\(ECM + \lambda{[\alpha\sum_{j}|\beta_j| +(1-\alpha)\sum_{j}\beta_j^2]}\)

\(\lambda\) controla toda la penalidad, mientras que \(\alpha\) controla la penalidad de la elastic net y tiende un puente entre Lasso y Ridge.

Con el parámetro \(\alpha\) indicamos con qué tipo de modelo deseamos trabajar:

Ridge: \(\alpha=0\)
Lasso: \(\alpha=1\)
Elastic Net: \(0<\alpha<1\)

Partición Train y Test

Realizamos una partición entre dataset de entrenamiento y testeo usando la función initial_split del paquete rsample.

train_test <- nba %>% initial_split(prop = 0.7)

train <- training(train_test)
test <- testing(train_test)

Lasso

En este caso vamos a trabajar con \(\alpha=1\).

¿Cuál es la penalización que introduce el modelo Lasso?
¿Cómo impacta esto en las variables?

# Vector con los salarios
nba_salary = train$salary
# Matriz con los regresores
nba_mtx = model.matrix(salary~., data = train)

# Modelo Lasso
lasso.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=1, # Indicador del tipo de regularizacion
                 standardize = F) # Que esta haciendo este parametro?
                 
lasso_coef = lasso.mod %>% tidy() %>% arrange(step)

lasso_coef

Gráficos de análisis

El comando plot nos permite realizar dos gráficos relevantes.

Gráfico de coeficientes en función del lambda

plot(lasso.mod, 'lambda')

Gráfico de coeficientes en función de la norma de penalización

plot(lasso.mod)

¿Qué muestra cada uno de estos gráficos?

Podemos realizar los gráficos para los valores de lambda también con ggplot.

g1=lasso_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
  geom_line() +
  theme_bw()  +
  theme(legend.position = 'none') +
  labs(title="Lasso con Intercepto",  y="Coeficientes")

g2=lasso_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
  geom_line() +
  theme_bw()  +
  theme(legend.position = 'none') +
  labs(title="Lasso sin Intercepto", y="Coeficientes")

plot_grid(g1,g2)

Veamos un poco mejor aquellas variables que sobreviven para mayores valores de lambda ¿Qué tienen en común todas estas variables?

# Seleccionamos los terminos que sobreviven para valores altos de lambda
terminos_sobrevientes = lasso_coef %>% 
  filter(log(lambda)>16.5, term != "(Intercept)") %>%
  select(term) %>% distinct() %>% pull()
# Graficamos
lasso_coef %>% filter(term %in% terminos_sobrevientes) %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
  geom_line(size=1) +
  geom_hline(yintercept = 0, linetype='dashed') +
  theme_bw() +
  labs(title="Lasso sin Intercepto", y="Coeficientes", subtitle= "\"Mejores\" variables") +
  scale_color_brewer(palette = 'Set1')

Vemos que las variables que “sobreviven” para mayores valores de lambda son las que están medidas con una escala mayor.

Estandarización en `glmnet`

Existen dos maneras de estandarizar las variables en glmnet.

Setear standardize = TRUE. Con esto se estandariza las regresoras y los coeficientes estimados estan en la escala original de la variable.
Pasar los conjuntos de datos estandarizados.

Lasso estandarizado

Replicamos todo el análisis previo pero para los datos estandarizados.

# Modelo lasso
lasso.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=1, # Indicador del tipo de regularizacion
                 standardize = TRUE) # Estandarizamos
                 
lasso_coef = lasso.mod %>% tidy() %>% arrange(step)

lasso_coef

Gráficos de análisis

Gráfico de coeficientes en función del lambda

plot(lasso.mod, 'lambda')

Gráfico de coeficientes en función de la norma de penalización

plot(lasso.mod)

Con ggplot

g1=lasso_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
  geom_line() +
  theme_bw()  +
  theme(legend.position = 'none') +
  labs(title="Lasso con Intercepto",  y="Coeficientes")

g2=lasso_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
  geom_line() +
  theme_bw()  +
  theme(legend.position = 'none') +
  labs(title="Lasso sin Intercepto", y="Coeficientes")

plot_grid(g1,g2)

Veamos ahora cuáles variables sobreviven para mayores valores de lambda.

# Seleccionamos los terminos que sobreviven para valores altos de lambda
terminos_sobrevientes = lasso_coef %>% filter(log(lambda)>13.1, term != "(Intercept)") %>%
  select(term) %>% distinct() %>% pull()
# Graficamos
lasso_coef %>% filter(term %in% terminos_sobrevientes) %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) +
  geom_line(size=1) +
  geom_hline(yintercept = 0, linetype='dashed') +
  theme_bw() +
  labs(title="Lasso sin Intercepto", y="Coeficientes", subtitle= "\"Mejores\" variables") +
  scale_color_brewer(palette = 'Set1')

Observamos que ahora tenemos otro set de “mejores” variables.

¿Podemos decidir cuál es el valor óptimo de lambda?

Elección del lambda óptimo

Para elegir el valor óptimo de lambda, lo común es realizar cross-validation. La función cv.glmnet nos permite realizar esto de manera sencilla.

Al igual que para la función glmnet cuenta con los parámetros:

x: matriz de variables
y: vector de la variable a predecir
alpha: tipo de modelo
standardize: flag lógico para estandarizar las variables

Nuevo parámetro:

type.measure: funció de pérdida/error que se va a utilizar en CV. Para los modelos de regularización el default es MSE.

Salida Base

lasso_cv = cv.glmnet(x=nba_mtx,y=nba_salary,alpha=1, standardize = T)
summary(lasso_cv)

##            Length Class  Mode     
## lambda     100    -none- numeric  
## cvm        100    -none- numeric  
## cvsd       100    -none- numeric  
## cvup       100    -none- numeric  
## cvlo       100    -none- numeric  
## nzero      100    -none- numeric  
## call         5    -none- call     
## name         1    -none- character
## glmnet.fit  12    elnet  list     
## lambda.min   1    -none- numeric  
## lambda.1se   1    -none- numeric

Brinda mucha información:

lambda: valor de lambda
cvm (Cross-validation mean): es la media del MSE (error)
cvsd (Cross-validation Standard Error): desvio estandar del MSE (error)
cvup y cvlo: Limite superior e inferior
nzero: Coeficientes distintos de cero
lambda.min: lambda para el cual el MSE (error) es mínimo
lambda.1se: lambda que se encuentra a 1 desvío estandar del lambda.min
glm.fit: incluye cantidad de variables, el valor de lambda y el porcentaje de deviance explicada por el modelo.

Si imprimimos el objeto tenemos:

lasso_cv

## 
## Call:  cv.glmnet(x = nba_mtx, y = nba_salary, alpha = 1, standardize = T) 
## 
## Measure: Mean-Squared Error 
## 
##      Lambda   Measure        SE Nonzero
## min  168429 3.626e+13 3.368e+12      17
## 1se 1188231 3.930e+13 3.814e+12       6

Gráfico Base

plot(lasso_cv)

El gráfico nos muestra la media del MSE con su límite superior e inferior y la cantidad de variables que sobreviven para cada valor de lambda.

Usando Broom

Obtenemos la información del objeto lasso_cv con las funciones tidy y glance.

# Información de CV en dataframe con tidy
lasso_cv %>% tidy()

# Lambda minimo y lambda a 1 desvio estandar
lasso_cv %>% glance()

Seleccionamos el lambda óptimo para crear el modelo final.

# Selección lambda óptimo
lasso_lambda_opt = lasso_cv$lambda.min

# Entrenamiento modelo óptimo
lasso_opt = glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=1, # Indicador del tipo de regularizacion
                 standardize = TRUE,  # Estandarizamos
                 lambda = lasso_lambda_opt)

# Salida estandar
lasso_opt

## 
## Call:  glmnet(x = nba_mtx, y = nba_salary, alpha = 1, lambda = lasso_lambda_opt,      standardize = TRUE) 
## 
##   Df %Dev Lambda
## 1 17 65.9 168400

# Tidy
lasso_opt %>% tidy()

# Glance (no es muy informativo)
lasso_opt %>% glance()

Han quedado 17 variables y el modelo explica el 65,9% de la deviance.

Ridge

En este caso vamos a trabajar con \(\alpha=0\). Vamos a replicar lo que ya realizamos para Lasso.

¿Cuál es la penalización que introduce el modelo Ridge?
¿Cómo impacta esto en las variables?

En este caso, ya seteamos parámetro para estandarizar las variables. Si no lo fijasemos, el default sería de igual modo standarize=TRUE.

#Modelo ridge
ridge.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0, # Indicador del tipo de regularizacion
                 standardize = TRUE)
#Coeficientes tidy                 
ridge_coef= ridge.mod %>% tidy()

ridge_coef

¿Qué ven de distinto en los coeficientes estimados del modelo respecto a Lasso?

Gráficos de análisis

Gráfico de coeficientes en función del lambda

plot(ridge.mod, 'lambda')

¿Qué ven de distinto en este gráfico respecto al que obtuvimos con la regresión Lasso?

Gráfico de coeficientes en función de la norma de penalización

plot(ridge.mod)

Gráficos de análisis con GGplot

g1=ridge_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Ridge con Intercepto",  y="Coeficientes")

g2=ridge_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Ridge sin Intercepto", y="Coeficientes")

plot_grid(g1,g2)

Elección del lambda óptimo

#cross-validation
ridge_cv=cv.glmnet(x=nba_mtx,y=nba_salary,alpha=0, standardize = T)

Gráfico Base

plot(ridge_cv)

Seleccionamos el lambda óptimo para crear el modelo final.

# Selección lambda óptimo
ridge_lambda_opt = ridge_cv$lambda.min

# Entrenamiento modelo óptimo
ridge_opt = glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0, # Indicador del tipo de regularizacion
                 standardize = TRUE,  # Estandarizamos
                 lambda = ridge_lambda_opt)

# Salida estandar
ridge_opt

## 
## Call:  glmnet(x = nba_mtx, y = nba_salary, alpha = 0, lambda = ridge_lambda_opt,      standardize = TRUE) 
## 
##   Df  %Dev  Lambda
## 1 50 65.91 2810000

# Tidy
ridge_opt %>% tidy() %>% mutate(estimate = round(estimate, 4))

En este caso conserva todas las variables, obteniendo un porcentaje de deviance explicado muy similar a Lasso: 65,9%.

Elastic Net

El modelo Elastic Net incorpora los dos tipos de penalización: Lasso (Norma L1) y Ridge (Norma L2). Define un compromiso entre las penalidades de lasso y ridge. El parámetro \(\alpha\) regula la importancia de cada penalización, cuanto más cerca de cero será más importante la penalización del tipo Ridge y más cerca de 1, la tipo Lasso.

En este caso vamos a trabajar con \(\alpha=0.5\). Vamos a replicar lo que ya realizamos para Lasso y Ridge, individualmente.

# Modelo elastic net
elastic.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0.5, # Indicador del tipo de regularizacion
                 standardize = TRUE)

# Coeficientes del modelo                 
elastic_coef= elastic.mod %>% tidy() %>% mutate(estimate = round(estimate, 4)) %>% arrange(step)  

elastic_coef

¿Qué ven de distinto en los coeficientes estimados del modelo respecto a Lasso y Ridge?

Gráficos de análisis

Gráfico de coeficientes en función del lambda

plot(elastic.mod, 'lambda')

¿Qué ven en este gráfico de distinto a los dos anteriores?

g1=elastic_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Elastic Net con Intercepto",  y="Coeficientes")

g2=elastic_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Elastic Net sin Intercepto", y="Coeficientes")

plot_grid(g1,g2)

Elección del lambda óptimo

elastic_cv=cv.glmnet(x=nba_mtx,y=nba_salary,alpha=0.5, standardize = T)

Grafico Base

Presten especial atención al eje superior ¿Qué está sucediendo?

plot(elastic_cv)

Seleccionamos el lambda óptimo para crear el modelo final.

# Selección lambda óptimo
elastic_lambda_opt = elastic_cv$lambda.min

# Entrenamiento modelo óptimo
elastic_opt = glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0.5, # Indicador del tipo de regularizacion
                 standardize = TRUE,  # Estandarizamos
                 lambda = elastic_lambda_opt)

# Salida estandar
elastic_opt

## 
## Call:  glmnet(x = nba_mtx, y = nba_salary, alpha = 0.5, lambda = elastic_lambda_opt,      standardize = TRUE) 
## 
##   Df  %Dev Lambda
## 1 19 65.37 405700

# Tidy
elastic_opt %>%  tidy()  %>% mutate(estimate = round(estimate, 4))

Breve comparación entre modelos

Vamos a comparar la relación entre el porcentaje de deviance explicada y lambda para los tres tipos de modelos que realizamos.

ridge_dev = ridge_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
  ggplot(., aes(log(lambda), dev.ratio)) +
  geom_point() +
  geom_line() +
  geom_vline(xintercept = log(ridge_lambda_opt), color='steelblue', size=1.5) +
  labs(title='Ridge: Deviance') +
  theme_bw() 

lasso_dev = lasso_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
  ggplot(., aes(log(lambda), dev.ratio)) +
  geom_point() +
  geom_line() +
  geom_vline(xintercept = log(lasso_lambda_opt), color='firebrick', size=1.5) +
  labs(title='Lasso: Deviance') +
  theme_bw()

elastic_dev = elastic_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
  ggplot(., aes(log(lambda), dev.ratio)) +
  geom_point() +
  geom_line() +
  geom_vline(xintercept = log(elastic_lambda_opt), color='forestgreen', size=1.5) +
  labs(title='Elastic Net: Deviance') +
  theme_bw()

plot_grid(ridge_dev, lasso_dev, elastic_dev)

Testing

Con los modelos óptimos que encontramos se puede calcular cuál es el RMSE en los datasets de training y testing, para decidir cuál es el modelo que minimiza el error en las predicciones.

# Definimos una función augment que funcione con glmnet
augment_glmnet = function(df,y, model) {
  # Matriz con los regresores
  formula = as.formula(str_c(y, "~.")) 
  data_matrix = model.matrix(formula, data = df)
  predictions = predict(model, data_matrix)
  pred_colname = str_c("predicted",y, sep="_")
  df[pred_colname] = predictions
  return(df)
}

Vamos a agregar dos modelos para comparar: el salario promedio del set de entrenamiento y el modelo lineal múltiple clásico.

# Salario promedio del set de entrenamiento
salario_promedio = mean(train$salary)
# Modelo lineal
modelo_lineal = lm(salary~., data = train)

Realizamos las predicciones ambos modelos.

# Salario promedio
prediccion_modelo_nulo = tibble(salary = train$salary, predicted_salary = salario_promedio)
# Predicciones del modelo lineal 
prediccion_modelo_lineal = augment(modelo_lineal) %>% mutate(predicted_salary = .fitted) %>% select(salary, predicted_salary)

Realizamos las predicciones de los modelos de glmnet.

# Lista de modelos
modelos_glmnet = list(lasso=lasso_opt, ridge=ridge_opt, elastic=elastic_opt)
# Predicciones de los modelos glmnet
lista_predicciones_training = map(.x = modelos_glmnet, .f = augment_glmnet, df=train, y="salary")
# Agregamos las otras predicciones a la lista 
lista_predicciones_training = lista_predicciones_training %>% prepend(list(nulo=prediccion_modelo_nulo, lineal = prediccion_modelo_lineal))

Obtenemos el RMSE para el set de entrenamiento.

map_dfr(.x = lista_predicciones_training, .f = rmse, truth="salary", estimate="predicted_salary", .id="modelo")

¿Cuál es el modelo que realiza la mejor predicción? ¿Qué esperan que suceda con el RMSE en el set de evaluación?

Obtenemos el RMSE para los modelos en el set de evaluación

# Predicción del promedio
prediccion_modelo_nulo = tibble(salary = test$salary, predicted_salary = salario_promedio)
# Predicción modelo lineal 
prediccion_modelo_lineal = augment(modelo_lineal, newdata = test) %>% mutate(predicted_salary = .fitted)
# Predicciones glmnet
lista_predicciones_test = map(.x = modelos_glmnet, .f = augment_glmnet, df=test, y="salary")
# Lista completa de predicciones
lista_predicciones_test = lista_predicciones_test %>% prepend(list(nulo=prediccion_modelo_nulo, lineal = prediccion_modelo_lineal))
# RMSE en el set de evaluación
map_dfr(lista_predicciones_test, rmse, truth="salary", estimate="predicted_salary", .id="modelo")

¿Cuál es el modelo de mejor performance en el set de evaluación? ¿Qué sucedió con el modelo de regresión clásico?

Regularización: Lasso, Ridge y Elastic Net

Juan Manuel Barriola y Sofia Perini

14 de Noviembre de 2020

Planteo del Problema

Conjuntos de datos

Análisis Exploratorio

Gráfico de la relación entre la posición y el salario

Correlograma

GGpairs (algunas variables)

Modelo Lineal

Coeficientes estimados

Evaluación de los modelos

Regularización

Partición Train y Test

Lasso

Gráficos de análisis

Estandarización en `glmnet`

Lasso estandarizado

Gráficos de análisis

Elección del lambda óptimo

Ridge

Gráficos de análisis

Gráficos de análisis con GGplot

Elección del lambda óptimo

Elastic Net

Gráficos de análisis

Elección del lambda óptimo

Breve comparación entre modelos

Testing

Regularización: Lasso, Ridge y Elastic Net

Juan Manuel Barriola y Sofia Perini

14 de Noviembre de 2020

Planteo del Problema

Conjuntos de datos

Análisis Exploratorio

Gráfico de la relación entre la posición y el salario

Correlograma

GGpairs (algunas variables)

Modelo Lineal

Coeficientes estimados

Evaluación de los modelos

Regularización

Partición Train y Test

Lasso

Gráficos de análisis

Estandarización en glmnet

Lasso estandarizado

Gráficos de análisis

Elección del lambda óptimo

Ridge

Gráficos de análisis

Gráficos de análisis con GGplot

Elección del lambda óptimo

Elastic Net

Gráficos de análisis

Elección del lambda óptimo

Breve comparación entre modelos

Testing

Estandarización en `glmnet`